Setup

Load libraries

library(tidyverse)
library(sf)
library(raster)
library(usdata)
library(tigris)
library(skimr)
library(DataExplorer)
library(raster)
library(terra)

Soil Erosion

RUSLE (A = R x K x LS x C x P)

Resolutions

  • Rainfall = 4 km
  • Erodibility = 30 m
  • Slope = 10 m
  • Crop scape = 30

Rainfall Erosivity

Load data

# rainfall - 2021
rnf <- stack("Data/Rainfall/dataset-sis-biodiversity-cmip5-global-83e8ca0e-4eea-4051-a911-c3299d0f2380 (1)/BIO12_gfdl-esm2m_rcp85_r1i1p1_1950-2100_v1.0.nc")

rnf <- mean(rnf[[72]])*3600*24*365*1000
plot(rnf)

Plot

#plot(rnf)
#sum_rnf <- sum(rnf)

plot(rnf)
plot(st_geometry(california), col = "red", add = TRUE)

Cropping and masking

Boundary lines

library(tigris)

# california
california <- counties(cb = TRUE, progress_bar = FALSE) %>% 
  filter(STATE_NAME == "California")

sj <- counties(cb = TRUE, progress_bar = FALSE) %>% 
  filter(NAME == "San Joaquin")

plot(st_geometry(california))
plot(st_geometry(sj), col = "blue", add = TRUE)

# cropping
crop_rnf <- raster::crop(rnf, extent(california))
plot(crop_rnf)

# masking
mask_rnf <- raster::mask(crop_rnf, mask = california)
plot(mask_rnf)

# resampling/interpolation
target <- mask_rnf
res(target) <- 0.0003176353
res_rnf <- raster::resample(mask_rnf, target, method = "bilinear")
plot(res_rnf)

raster::writeRaster(res_rnf, filename = file.path("Data/Rainfall/PRISM//res_rnf_california.tif"), format="GTiff", overwrite=TRUE)
res_rnf <- raster("Data/Rainfall/PRISM/res_rnf_california.tif")
plot(res_rnf)
plot(st_geometry(sj), add = TRUE)

# cropping
crop_rnf <- raster::crop(res_rnf, extent(sj))
plot(crop_rnf)

# masking
mask_rnf <- raster::mask(crop_rnf, mask = sj)
plot(mask_rnf)

raster::writeRaster(mask_rnf, filename = file.path("Data/Rainfall/PRISM//res_rnf_sj.tif"), format="GTiff", overwrite=TRUE)

Erosivity factor

# load data 
rnf_sj <- raster("Data/Rainfall/res_rnf_sj.tif")
plot(rnf_sj)

# The equation of Moore
# Implement the equation of Moore as a function in R
calculate_r_moore <- function(p) {
  ke <- 11.46*p - 2226
  r <- 0.029*ke - 26
  r_si <- 17.02*r # Conversion from imperial to SI units
  return(r_si)
}

r_moore <- calc(x = rnf_sj, fun = calculate_r_moore)
plot(r_moore)

# save raster
raster::writeRaster(r_moore, filename = file.path("Data/Rainfall/erosivity_factor.tif"), format="GTiff", overwrite=TRUE)

Erodibility factor (K)

Load data

ef <- raster("Data/Erodibility/SJ_Erodibility.tif")
raster::plot(ef)

Slope length and steepness

DEM

Load data

dem1 <- raster("Data/DEM or DSM/USGS_13_n38w121_20220103.tif")
dem2 <- raster("Data/DEM or DSM/USGS_13_n38w122_20220810.tif")
dem3 <- raster("Data/DEM or DSM/USGS_13_n39w121_20220206.tif")
dem4 <- raster("Data/DEM or DSM/USGS_13_n39w122_20220206.tif")

dem_2022 <- raster::mosaic(dem1, dem2, dem3, dem4,  fun = 'mean')
plot(dem_2022)
#plot(st_geometry(california), add = TRUE)
plot(st_geometry(sj), add = TRUE, col = "red")


sj =sf::st_set_crs(sj, crs(dem_2022))
plot(dem)
plot(st_geometry(sj), add = TRUE)

Mosaic

dem_2022 <- raster::mosaic(dem1, dem2, dem3, dem4,  fun = 'mean')

# plot DEM
plot(dem_2022)
plot(st_geometry(sj), add = TRUE, col = "red")

Crop and Mask

# crop
crop_dem <- crop(dem_2022, extent(sj))
plot(crop_dem)

# mask
mask_dem <- mask(crop_dem, sj)
plot(mask_dem)

# save 
raster::writeRaster(mask_dem, filename = file.path("Data/DEM or DSM/sj_dem.tif"), format="GTiff", overwrite=TRUE)

Load clipped dem

sj_dem <- raster("Data/DEM or DSM/sj_dem.tif")
plot(sj_dem)

Slope

slope <- terrain(sj_dem, "slope", unit = "degrees")
plot(slope)

# save 
#raster::writeRaster(slope, filename = file.path("Data/DEM or DSM/sj_slope.tif"), format="GTiff", overwrite=TRUE)

LS from ArcGIS Pro

# flow accumulation
ls <- raster("Data/DEM or DSM/LS_1_1/LS_1_1.tif")

plot(ls)
plot(st_geometry(sj), add = TRUE)

# crop
crop_ls <- crop(ls, extent(sj))
mask_ls <- mask(crop_ls, sj)
plot(mask_ls)

mask_ls
ef

# resampling/interpolation
target <- slope
res(target) <- res(ef) # 30 meters
mask_ls <- raster::resample(mask_ls, target, method = "bilinear")
plot(mask_ls)

# save 
raster::writeRaster(mask_ls, filename = file.path("Data/LS /LS_factor.tif"), format="GTiff", overwrite=TRUE)
# load LS factor
ls_factor <- raster("Data/LS Factor/LS_factor.tif")
plot(ls_factor)

Cover management

Load data

cs <- raster("Data/Crop cover/Reclassify_SJ_Cropcover_2021.tif")
plot(cs)

# crop and mask
crop_cs <- crop(cs, extent(sj))
mask_cs <- mask(crop_cs, sj)
plot(mask_cs)

# crop factor
c_factor <- mask_cs/100
plot(c_factor)

Support Practices: Cropping (P-factor)

  • Cmanagementpractices = 1 for conventional Strip cropping;
  • Cmanagementpractices = 0.35 for conservation contour cropping;
  • Cmanagementpractices = 0.25 for terrace cropping.

Erosion calculation

setup

e_factor <- raster::resample(ef, r_moore, method = "bilinear")
ls_factor <- raster::resample(ls_factor, r_moore, method = "bilinear")
c_factor <- raster::resample(c_factor, r_moore, method = "bilinear")

First scenario

Conventional Strip Cropping

RUSLE_1 <- r_moore * e_factor * ls_factor * c_factor * 1
plot(RUSLE_1)

# save 
raster::writeRaster(RUSLE_1, filename = file.path("Data/RUSLE/RUSLE_1.tif"), format="GTiff", overwrite=TRUE)

Second scenario

Conservation Contour Cropping

RUSLE_2 <- r_moore * e_factor * ls_factor * c_factor * 0.35
plot(RUSLE_2)

# save 
raster::writeRaster(RUSLE_2, filename = file.path("Data/RUSLE/RUSLE_2.tif"), format="GTiff", overwrite=TRUE)

Third scenario

Terrace Cropping

RUSLE_3 <- r_moore * e_factor * ls_factor * c_factor * 0.25
plot(RUSLE_3)

# save 
raster::writeRaster(RUSLE_3, filename = file.path("Data/RUSLE/RUSLE_3.tif"), format="GTiff", overwrite=TRUE)

Comparison of scenarios

RUSLE 1 vs RUSLE 2

r1_r2 <- RUSLE_1 - RUSLE_2
plot(r1_r2)

Projections

library(wesanderson)
# Gradient color
pal <- wes_palette("Zissou1", 100, type = "continuous")

Load ppt data

fut_ppt <- stack("Data/Rainfall/dataset-sis-biodiversity-cmip5-global-83e8ca0e-4eea-4051-a911-c3299d0f2380 (1)/BIO12_gfdl-esm2m_rcp85_r1i1p1_1950-2100_v1.0.nc")
plot(fut_ppt, col = pal)

Referencing

library(tigris)
us <- states(cb = TRUE) %>% 
  filter(!NAME %in% c("Puerto Rico", "Guam", "American Samoa", "United States Virgin Islands", "Commonwealth of the Northern Mariana Islands", "Alaska", "Hawaii")) 
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |======================================================================| 100%
#plot(fut_ppt[[1]])
plot(st_geometry(us))

Select study years

  • 2030

  • 2050

  • 2070

  • 2100

# 2030s
y2030 <- fut_ppt[[81]]*3600*24*365*1000
plot(y2030)

# 2050s
y2050 <- fut_ppt[[101]]*3600*24*365*1000
plot(y2050)

# 2070s
y2070 <- fut_ppt[[121]]*3600*24*365*1000
plot(y2070)

# 2100s
y2100 <- fut_ppt[[151]]*3600*24*365*1000
plot(y2100)

# stack
fut_ppt <- stack(y2030, y2050, y2070, y2100)
names(fut_ppt) <- c("2030", "2050", "2070", "2100")
plot(fut_ppt)

Cropping and Masking

# cropping
crop_ppt <- raster::crop(fut_ppt, extent(california))
plot(crop_ppt)

# masking
mask_ppt <- raster::mask(crop_ppt, mask = california)
plot(mask_ppt)

# load target layer
res_rnf <- raster("Data/Rainfall/PRISM/res_rnf_california.tif")

# resampling/interpolation
res_ppt <- raster::resample(mask_ppt, res_rnf, method = "bilinear")
plot(res_ppt)

raster::writeRaster(res_ppt, filename = file.path("Data/FUTURE_PROJECTIONS/res_ppt_california.tif"), format="GTiff", overwrite=TRUE)
res_ppt <- brick("Data/FUTURE_PROJECTIONS/res_ppt_california.tif")
plot(res_ppt)

plot(res_ppt[[1]])
plot(st_geometry(sj), add = TRUE)

# cropping
crop_ppt <- raster::crop(res_ppt, extent(sj))
plot(crop_ppt)

# masking
mask_ppt <- raster::mask(crop_ppt, mask = sj)
plot(mask_ppt)

raster::writeRaster(mask_ppt, filename = file.path("Data/FUTURE_PROJECTIONS/res_ppt_sj.tif"), format="GTiff", overwrite=TRUE)
res_ppt <- stack("Data/Future projections/res_ppt_sj.tif")
#res_ppt <- res_ppt*3600*24*365*1000
names(res_ppt) <- c("2030", "2050", "2070", "2100")
plot(res_ppt)

# comparison of decades
ppt_df <- res_ppt %>% 
  as.data.frame(xy = TRUE) %>% 
  pivot_longer(names_to = "key", values_to = "value", cols = c(3:6)) %>% 
  filter(!is.na(value))

ppt_df %>% 
  ggplot() +
  geom_raster(aes(x = x, y = y, fill = value)) +
  scale_fill_distiller(palette = "RdPu", direction = 1) +
  theme_minimal() +
  facet_wrap(~ key)

# erosivity factor
r_factor_future <- calc(x = res_ppt, fun = calculate_r_moore)

# setup
r_factor_future <- raster::resample(r_factor_future, r_moore, method = "bilinear")
plot(r_factor_future)

First scenario

Conventional Strip Cropping

RUSLE_1_future <- r_factor_future * e_factor * ls_factor * c_factor * 1
names(RUSLE_1_future) <- c("2030", "2050", "2070", "2100")
plot(RUSLE_1_future)

# save 
raster::writeRaster(RUSLE_1_future, filename = file.path("Data/RUSLE/RUSLE_1_future.tif"), format="GTiff", bylayer = TRUE, overwrite=TRUE)

Second scenario

Conservation Contour Cropping

RUSLE_2_future <- r_factor_future * e_factor * ls_factor * c_factor * 0.35
names(RUSLE_2_future) <- c("2030", "2050", "2070", "2100")
plot(RUSLE_2_future)

# save 
raster::writeRaster(RUSLE_2_future, filename = file.path("Data/RUSLE/RUSLE_2_future.tif"), format="GTiff", bylayer = TRUE, overwrite=TRUE)

Third scenario

Terrace Cropping

RUSLE_3_future <- r_factor_future * e_factor * ls_factor * c_factor * 0.25
names(RUSLE_3_future) <- c("2030", "2050", "2070", "2100")
plot(RUSLE_3_future)

# save 
raster::writeRaster(RUSLE_3_future, filename = file.path("Data/RUSLE/RUSLE_3_future.tif"), format="GTiff", bylayer = TRUE, overwrite=TRUE)

EDA

Rainfall

rnf_sj <- raster::resample(rnf_sj, res_ppt, method = "bilinear")
names(rnf_sj) <- "2021"

plot(stack(rnf_sj, res_ppt), col = rev(pal))